BIO4000W GIS Deliverable

View the html here: https://raw.githack.com/jess-devine/BIO4000W_GIS/refs/heads/main/README.html

Is Southern Afrotemperate Forest Expanding in the Cape Peninsula?

Investigating whether pioneer forest species indicates forest expansion. Using GIS and species occurrence data to assess forest dynamics.

Species that are Southern Afrotempraet Forest pioneers: Virgilia oroboides (Cape Keurboom), Virgilia divaricata (Garden Route Keurboom), Rapanea melanophloeos (Cape Beech), and Kiggelaria africana (Wild Peach).

Observation data from iNaturalist and City of Cape Town vegetation layers will be used.

## Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
## 
## Attaching package: 'gridExtra'
## 
## 
## The following object is masked from 'package:dplyr':
## 
##     combine
## 
## 
## 
## Attaching package: 'cowplot'
## 
## 
## The following object is masked from 'package:patchwork':
## 
##     align_plots
## 
## 
## The following object is masked from 'package:lubridate':
## 
##     stamp
## 
## 
## Loading required package: viridisLite
species_list <- c("Virgilia oroboides", "Virgilia divaricata", 
                  "Rapanea melanophloeos", "Kiggelaria africana")
species_plots <- list()
species_data <- list()

# get veg data 
veg <- st_read("data/cape_peninsula/cape_peninsula/veg/Vegetation_Indigenous.shp")
## Reading layer `Vegetation_Indigenous' from data source 
##   `/home/jess/GIT/BIO4000W_GIS/data/cape_peninsula/cape_peninsula/veg/Vegetation_Indigenous.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1325 features and 5 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -63972.95 ymin: -3803535 xmax: 430.8125 ymax: -3705149
## Projected CRS: WGS_1984_Transverse_Mercator
st_write(veg, "data/cape_peninsula/cape_peninsula/veg/Vegetation_Indigenous_duplicate.shp", append = FALSE)
## Deleting layer `Vegetation_Indigenous_duplicate' using driver `ESRI Shapefile'
## Writing layer `Vegetation_Indigenous_duplicate' to data source 
##   `data/cape_peninsula/cape_peninsula/veg/Vegetation_Indigenous_duplicate.shp' using driver `ESRI Shapefile'
## Writing 1325 features with 5 fields and geometry type Polygon.
vegr <- st_read("data/cape_peninsula/cape_peninsula/veg/Vegetation_Indigenous_Remnants.shp")
## Reading layer `Vegetation_Indigenous_Remnants' from data source 
##   `/home/jess/GIT/BIO4000W_GIS/data/cape_peninsula/cape_peninsula/veg/Vegetation_Indigenous_Remnants.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 3428 features and 7 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -63951.23 ymin: -3803532 xmax: 420.7595 ymax: -3705506
## Projected CRS: WGS_1984_Transverse_Mercator
st_write(vegr, "data/cape_peninsula/cape_peninsula/veg/Vegetation_Indigenous_duplicate.shp", append = FALSE)
## Deleting layer `Vegetation_Indigenous_duplicate' using driver `ESRI Shapefile'
## Writing layer `Vegetation_Indigenous_duplicate' to data source 
##   `data/cape_peninsula/cape_peninsula/veg/Vegetation_Indigenous_duplicate.shp' using driver `ESRI Shapefile'
## Writing 3428 features with 7 fields and geometry type Polygon.
# Crop vegr 
  ext <- c(-66642.18, -3809853.29, -44412.18, -3750723.29)
  names(ext) <- c("xmin", "ymin", "xmax", "ymax") 
  vegr <- st_crop(vegr, ext)
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
for (species_name in species_list) {
  # Retrieve iNaturalist data
  inat_data <- get_inat_obs(taxon_name = species_name,
                            bounds = c(-35, 18, -33.5, 18.5),
                            maxresults = 1000)
  
  # Filter iNaturalist data
  inat_data <- inat_data %>% 
    filter(positional_accuracy < 46 & 
           latitude < 0 &
           !is.na(latitude) &
           captive_cultivated == "false" &
           quality_grade == "research")

  # Convert to spatial object
  inat_data <- st_as_sf(inat_data, coords = c("longitude", "latitude"), crs = 4326)
  inat_data <- st_transform(inat_data, st_crs(vegr))
  inat_data <- st_intersection(inat_data, vegr)  # Remove points in urban areas

  # Ensure National_ is a character (fix color mapping issue)
  inat_data$National_ <- as.character(inat_data$National_)
  
  # Store in the list
  species_data[[species_name]] <- inat_data
  
  # buffer points not in Southern Afrotemprate Forest
  nsaf <- inat_data %>% 
  filter(National_ != "Southern Afrotemperate Forest") %>%
  st_buffer(dist = 250) 
  #Intersect new polygons with veg remnants and filter for those that overlap Southern Afrotemperate Forest only
  nsaf <- st_intersection(nsaf, vegr) %>% filter(National_.1 == "Southern Afrotemperate Forest")

  # Get count of non-buffered points per vegetation type
  no_buffer_summary <- inat_data %>%
    st_drop_geometry() %>%
    group_by(National_) %>%
    summarise(count_no_buffer = n())

  # Get count of buffered points per vegetation type
  buffer_change_summary <- nsaf %>%
    st_drop_geometry() %>%
    group_by(National_) %>%
    summarise(count_buffer_change = n())  # Renamed to avoid confusion

  # Compute total buffered points
  total_buffered <- sum(buffer_change_summary$count_buffer_change, na.rm = TRUE)

  # Merge summaries
  combined_summary <- full_join(no_buffer_summary, buffer_change_summary, by = "National_")
  combined_summary <- combined_summary %>%
    mutate(count_no_buffer = replace_na(count_no_buffer, 0),
           count_buffer_change = replace_na(count_buffer_change, 0))

  # Adjust count_buffer correctly
  combined_summary <- combined_summary %>%
    mutate(count_buffer = ifelse(National_ == "Southern Afrotemperate Forest",
                                 count_no_buffer + total_buffered,  # SAF gets the initial count plus total buffer count 
                                 pmax(0,count_no_buffer - count_buffer_change)))  # Others get adjusted values
  
  # Rename columns for clarity
combined_summary <- combined_summary %>%
  rename(
    "Vegetation Type" = National_,
    "Initial Count" = count_no_buffer,
    "Buffered Change" = count_buffer_change,
    "Final Count" = count_buffer
  )

  # Create the Plot
  species_plot <- ggplot() +
    annotation_map_tile(type = "osm", progress = "none") +
    geom_sf(data = inat_data, aes(color = National_), shape = 16) +  # Uniform shape
    labs(title = paste("Observations of", species_name)) +
    theme_minimal()

  # Create the Table
  table_plot <- tableGrob(combined_summary)

  # Print Separately
  print(species_plot)  # Print the figure

  # Move to a new page & draw the table
  grid.newpage()
  grid.draw(table_plot)  
}
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries

## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries

## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries

## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries

All of the tree species, except Rapanea melanphloeos, had more observations more than 250m from Southern Afrotemprate Forest than within this vegetation type. However, the extent of Southern Afrotemprate Forest is small on the Cape Peninsula.

The majority of species seem to be expanding into Peninsula Sandstone Fynbos more than other vegetation types. However, Kigillaria africana seems to be expanding more into Peninsula Granite Fynbos.

# Combine all species data into a single dataframe
combined_data <- bind_rows(species_data, .id = "Species")

# Define a color palette for vegetation type
veg_colors <- viridis::viridis(length(unique(combined_data$National_)))

# Define unique shapes for each species
species_shapes <- c("Virgilia oroboides" = 16,  # Filled circle
                    "Virgilia divaricata" = 17, # Filled triangle
                    "Rapanea melanophloeos" = 18, # Filled diamond
                    "Kiggelaria africana" = 15) # Filled square

# Generate the final plot
ggplot() +
  annotation_map_tile(type = "osm", progress = "none") +
  geom_sf(data = combined_data, 
          aes(color = National_, shape = Species), 
          size = 3) +  # Adjust size for better visibility
  scale_color_manual(values = veg_colors) +  # Assign vegetation colors
  scale_shape_manual(values = species_shapes) +  # Assign species shapes
  labs(title = "All Four Species Observations", 
       color = "National Vegetation Type", 
       shape = "Species") +
  theme_minimal() +
  theme(legend.position = "right")

# Ensure data is in WGS84 (longitude-latitude)
combined_data <- combined_data %>% st_transform(4326)
vegr <- st_transform(vegr, 4326)

# Define a color palette for speceis
library(RColorBrewer)
species_palette <- colorFactor(brewer.pal(min(9, length(unique(combined_data$Species))), "Set1"), combined_data$Species)

# Define a color palette for vegetation type
veg_palette <- colorFactor(viridis(length(unique(vegr$National_))), vegr$National_)

# Create iNaturalist links for popups
combined_data <- combined_data %>%
  mutate(iNat_link = paste0("<a href='https://www.inaturalist.org/observations/", id, 
                            "' target='_blank'>View on iNaturalist</a>"))

# Leaflet Map
leaflet(combined_data) %>%
  addTiles() %>%
  addPolygons(
    data = vegr,
    fillColor = ~veg_palette(National_),  # Use colorFactor function
    color = ~veg_palette(National_), 
    weight = 1,
    fillOpacity = 0.5,
    popup = ~paste("<b>Vegetation Type:</b>", National_),
    group = "Vegetation Remnants"
  ) %>%
  addCircleMarkers(
    data = combined_data,
    lng = ~st_coordinates(geometry)[,1], 
    lat = ~st_coordinates(geometry)[,2], 
    radius = 5, 
    color = ~species_palette(Species),  # Use colorFactor function
    fillOpacity = 0.9,
    popup = ~paste("<b>Species:</b>", Species, "<br>",
                   "<b>Vegetation Type:</b>", National_, "<br>",
                   iNat_link),
    group = "Observations"
  ) %>%
  addLegend("bottomright", 
            pal = species_palette, 
            values = combined_data$Species, 
            title = "Species") %>%
  addLegend("topright", 
            pal = veg_palette, 
            values = vegr$National_, 
            title = "Vegetation Type") %>%
  addLayersControl(
    overlayGroups = c("Vegetation Remnants", "Observations"),
    options = layersControlOptions(collapsed = FALSE)
  )

From these maps it apears that Southern Afrotemprate forest may be expanding in the Northern half of the Cape Peninsula. Historical pine plantations that are now nature reserves in this region may be hotspots for forest expansion, likely because fire was excluded for a long period of time while the plantations were in operation. Additionally, Silvermine dam may be facilitating Southern Afrotemprate Forest expansion.

However, some of the observations are in parks, botanical gardens and plantations. This likely accounts for the majority of points in Cape Flats Sand Fynbos. Aditionally, it is likely natural for there to be Southern Afrotemprate Forest to exit around small streams and rivers but the vegetation map scal is too course to map these.